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Abstract 

The behaviour of a solid-liquid-gas system near the three-phase contact line is consid- 
ered using a diffuse-interface model with no-slip at the solid and where the fluid phase 
is specified by a continuous density field. Relaxation of the classical approach of a sharp 
liquid-gas interface and careful examination of the asymptotic behaviour as the contact 
line is approached is shown to resolve the stress and pressure singularities associated with 
the moving contact line problem. Various features of the model are scrutinised, along- 
side extensions to incorporate slip, finite-time relaxation of the chemical potential, or a 
precursor film at the wall. 

1 Introduction 

A moving contact line occurs at the location where two ostensibly immiscible fluids and a 
solid meet. It arises in a wide range of both natural and technological processes, from insects 
walking on water [15j and the wetting properties of plant leaves jT3] to coating [51], inkjet 
printing [8, 38 and oil recovery [27j. In addition to its crucial role in wide-ranging applications, 
it remains a persistent problem, a long-standing and fundamental challenge in the field of fluid 
dynamics, despite its apparent simplicity at first sight (see e.g. review articles [TOl \7\ l2j [3]). 
Not surprisingly it has been investigated extensively, both experimentally and theoretically, 
for several decades. 

At the heart of the moving contact line problem is that, when treated classically as 
two immiscible fluids moving along a solid surface satisfying the no-slip condition, there is 
no solution due to the multivalued velocity at the contact line, [11, 33J. This is known 
most famously in the literature as the problem of a non-integrable stress singularity, a result 
published a few decades ago along with the nonphysical prediction that an infinite force is 
required to submerge a solid object Hi*. 

The resolution of the problem may have initially appeared straightforward. The no-slip 
condition at the wall could not be satisfied, thus some form of slip in the contact line vicinity 
should be allowed. Navier-slip, written down in the early 19th Century [19], was a prime 
candidate, a form of which was suggested in the concluding remarks of [16]. The fact that 
wetting and the moving contact line remain an open debate and a fruitful research area is 



largely due to the particular microscale ingredients that may alleviate the problem being nu- 
merous and hotly debated — see e.g. the wide range of discussion articles recently, |36j . Various 
alternative models to slip at the wall include: a precursor film ahead of the contact line [30J; 
rheological effects |4l)j ; treating surfaces as separate thermodynamic entities with dynamic 
surface tensions [32J (see a recent critical investigation of this model, |35|): introducing nu- 
merical slip [26]; including evaporative fluxes [6]; and considering the interface to be diffuse, 
numerical work reviewed e.g. in pQ. 

In this Letter, we examine analytically a diffuse-interface model without any other in- 
gredients, being both self-consistent and physically relevant: rather than considering a sharp 
fluid-fluid interface as a surface of zero thickness where quantities {e.g. density) are, in gen- 
eral, discontinuous, considers the interface to have a non-zero thickness with quantities varying 
smoothly but rapidly, in agreement with developments from the statistical mechanics of liquids 
{e.g. \\.2\ 144] ) which show that interfaces have a finite width of the order of several molecular 
diameters. The fluid density p then acts as an order parameter such that in the sharp-interface 
limit the two bulk fluids satisfy p = pi and p = py = 0, being liquid and vapour respectively, 
where we consider the behaviour of the system with vapour phase of negligible density. The 
liquid-gas interface is then defined as the location where p = {pl + Pv)/2- 

Diffuse-interface models have been popular for numerical implementation as tracking of 
the fluid-fluid interface is not required in the resulting free-boundary problem, instead the 
interface is inferred from density field contours. For solid- liquid- gas systems the seminal study 
of Seppecher [31] is often referred to when suggesting that diffuse-interface models resolve 
the moving contact line problem. Whilst this contains some discussion of the asymptotics, 
the analysis was largely incomplete, with asymptotic regions being probed without careful 
justification and the crucial behaviour close to the contact line only investigated numerically (a 
number of constraints were also imposed, e.g. 90° contact angles and fluids of equal viscosity). 
Full numerical simulations for the liquid-gas problem have also be undertaken in other studies 
{e.g. [H [5]), and binary fluids by a number of authors using diffuse-interface methods of a 
different form, where a coupled Cahn-Hilliard equation models the diffusion between the two 
components [TH1 19| ITT] . 

Here, we undertake an analytical investigation by considering the contact line behaviour 
for a liquid-gas system with two basic elements: (a) the interface has a finite thickness, which 
is expected from statistical mechanics studies as noted earlier, and (b) the no-slip condition is 
applied at the wall. This diffuse-interface model then resolves the moving contact line problem 
without the need for any additional physics. An important ramification of this analysis is 
that the wetting boundary condition used in conjunction with diffuse-interface in existing 
numerical studies needs to be appropriately modified, otherwise it leads to a density gradient 
near the wall at large distances from the contact line. The possibility of density variations 
such as these are often included when disjoining pressure models are considered, where many 
studies utilise the long- wave (or lubrication) approximation, e.g. [71 [20]. A recent study of 
intermolecular forces in the contact line region using approaches from statistical mechanics, 
namely density- functional theory [21], demonstrated that for partially wetting fluids there 
is no precursor film on the microscale, but the theory predicts a nanoscale film of a few 
molecular diameters adsorbed in front of the contact line. In the macro/mesoscopic setting 
of the current work, we thus predominately consider the case where the bulk densities are 
valid up to the walls (but other cases are considered), and we make no assumptions on thin 
films — rather considering arbitrary finite contact angles. 
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Figure 1: Sketch of the diffuse- interface model near a wall. 

2 Problem specification 

Consider a fluid in the upper half (x, y)-plane Q with a solid surface dQ at y = 0; see Fig. [TJ 
The free energy of the system has contributions from an isothermal fluid with a Helmholtz free 
energy functional and from a wall energy f w = f w (p), thus given by = £df2 + f dn f w dA, 
where 

£ = pf{p) + K\Vp\ 2 /2 - Gp, (1) 

with G the chemical potential, K a gradient energy coefficient (assumed constant for simplic- 
ity), and pf(p) a double well potential chosen to give the two equilibrium states p = {O,^}. 
Such a form for the free energy and associated diffuse-interface approximations has been used 
by numerous authors for wetting problems such as Seppecher [31], Pismen and Pomeau [23], 
and Pismen [22] : see also the review by de Gennes pTJ. The effect of the non-local terms 
neglected in the local approximation to obtain such a free energy has been considered at 
equilibrium in the aforementioned studies [22\ [2~T] , where the long-range intermolecular inter- 
actions are responsible for an algebraic decay of the density profile away from the interface 
instead of the exponential one as predicted here (seen later, in eq. ([7])). In our dynamic 
situation, we focus on the local approximation to elucidate the contact line behaviour in 
a simplified, yet widely used setting. The density field augments the usual hydrodynamic 
equations via the capillary (or Korteweg) stress tensor T through 

where I is the identity tensor, and with e being the interface width, and T arising from 
Noether's theorem, pQ. Using eq. ([1]) and coupling in the compressible Stokes equations 
(assuming creeping flow) the capillary tensor with the usual viscous stress tensor r, taken as 
deviatoric for simplicity, yields 

dtp + V-(pu) = 0, V-(T + r) = 0, 
T = (pf(p) + K\ Vp| 2 /2 - Gp) I - KVp Vp, 
t = p(p) [(Vu) + (Vu) T - 2(V • u)I/3] , 

G = -KV 2 p + dp(pf(p)), 

where u and p(p) are the fluid velocity and viscosity, respectively, t is time, and the thermo- 
dynamic pressure is given by p = p 2 f'(p). The form for G arises from the Euler-Lagrange 
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equation corresponding to the free energy, and we take p(p) 
for the two equilibrium states as pi and psy = 0. 

On the solid surface d£l, with normal n„,, we impose 



= hlp/ PL, giving the viscosities 



u = u 



■w ! 



Kn w -Vp + fUp) = 0, 



with wall velocity = (— V, 0) in Cartesian coordinates. The first condition is the classical 
no-slip, whilst the second arises through variational arguments and is termed the natural 
(or wetting) boundary condition. The form for f w should be chosen to satisfy Young's law 
at the contact line, with solid-liquid, solid-gas and liquid- gas surface tensions f w (pL) = 
fw(pv) = &v and o respectively, and with contact angle 9s- A cubic is the lowest-order 
polynomial required such that the wall free energy can be minimised for the bulk densities 
and prevent a film forming ahead of the contact line, i.e. f' w {p£) = f' w {fy = 0- This is unlike 
the linear forms proposed in the previous studies for liquid-gas problems [3TJ HI EJ 03] , allowing 
us to consider a diffuse-interface model without any additional physics from the microscale 
(although this will be relaxed in the following section). We define 



fw{p) = [p L Vcos^ 5 (4p 3 - 6p 2 p L + Pl) + <?v + °l\ A 
giving f w (p) = —Qap(pL — p) cos 9 s/ p\, and f w (0) - J w {pl) = crcos9 s , with Young's law 



thus satisfied. It is noteworthy that the natural boundary condition may be replaced by a 
constant density condition if a precursor film/ disjoining pressure model is to be considered, 
i.e. p = p a on dfl (as used in [23J, and considered in the following section), although as we 
shall demonstrate, this model is not necessary. Finally, for a one-dimensional density profile 
p{z) in equilibrium ((5 = 0, u = 0), the surface tension across the interface is 



To nondimensionalise, we use typical length, velocity and density scales X, V and pi 
respectively. The pressure and viscous stress are scaled with p^V/X, and the capillary stress 
with Kp 2 L /(eX). Finally, G is scaled with Kp L /{eX) and / with Kp L /e 2 . The governing 
equations then contain the nondimensional parameters 



being the Cahn number and a modified Capillary number based on the model parameter 
K, respectively. Nondimensional variables are denoted as their dimensional counterparts 
with bars dropped, e.g. nondimensional bulk densities are p = {0, 1}. Ca^ is related to the 
usual Capillary number, Ca = p£,V/a = 6Ca/c, through eq. ([2]). The governing equations in 
nondimensional form are 




(2) 



Cn = e/X, and Ca^ = pLVe/{Kp L ) 



d t p + V • (pu) = 0, M = Cafc T + r, V • M = 0, 
T = (Cn _1 /9/(p) + Cn| Vp| 2 /2 — Gp) I - Cn Vp ® Vp, 
r = p [(Vu) + (Vu) T - 2(V • u)I/3] , 

G = -CnV 2 p + Cn- 1 d p (pf(p)), 



(3) 



where M is introduced as the total stress tensor, p = (CnCa^) 1 p 2 f'(p), and 



p / \ p M x2 ■ . P 2 (l - 3/C 

/(P) = |(l-P), givmg P= 2C ^ 




(4) 
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Figure 2: The equilibrium density behaviour (contours) near the contact line for 8 s = tt/4, 
from eq. ([7]) in inner variables. 

On the solid surface <9f2, we have 

u = u w , Cnn„ • Vp = cos6» s (l - p)p, (5) 

where u w = (—1,0), and n w = (0, —1), in Cartesian components. We can rewrite the govern- 
ing equations as 

Ca^ 1 Cn pV(V 2 p) - Vp + p[V 2 u + V(V • u)/3] 

+ [(Vu) + (Vu) T - 2(V • u)I/3]Vp = 0, (6) 

where Vp = (2 Cn Ca^) -1 V[p 2 (l — 3p)(l — p)], with the (steady) momentum equation V • 
(pu) = 0, and the boundary conditions on d£l remaining as above, in eq. ([5]). We initially 
consider the equilibrium behaviour of the system, corresponding to Ca^ <C 1, where eqs. (0)- 
are thus reduced to 

2Cn 2 pV(V 2 p) = V [p 2 (l - 3p)(l - p)] , 

subject to — Cndyp = cos #s(l — p)p at y = 0, and with p — > {0, 1} and V 2 p — > as x — > ±oo 
to obtain the expected bulk behaviour. The solution subject to the above conditions is 

p= (1 -tanh [(2Cn)- 1 (xsin6»5 + ycos6»5)]) /2, (7) 

having also fixed the interface at p = 1/2. This profile is flat and at angle 9s, shown in Fig. [2] 
in inner variables (where {x,y} = Cn{z,y}, for comparison to forthcoming plots). 

For physical systems, the scale over which density varies between liquid and gas is much 
smaller than the macroscopic length scale, and hence Cn 1. The asymptotic behaviour as 
Cn — >• is known as the sharp-interface limit, and understanding of it is of vital importance 
when considering diffuse-interface models as classical continuum models should be recovered 
if correct predictions are to be found. 

A careful asymptotic analysis of the outer solution away from the interface, and of the 
interfacial region away from the wall (using body fitted coordinates), shows that the expected 
sharp-interface equations (the Stokes equations with no-slip and the usual capillary surface 
stress conditions) are indeed recovered. A recent study for general phase-field models of 
interfaces without walls shows that the classical model is only found in the sharp-interface 
limit when at equilibrium, |13] , This is not the case here. 




5 



Consider now the inner region near to both interface and wall in polar coordinates with 
r = O(Cn). The scaling r = Cnf (inner variables denoted with tildes) retains all terms in 
the governing equations and boundary conditions, giving a complete dominant balance. Of 
particular interest is the behaviour as the contact line is approached — the location where a 
stress singularity or no solution (due to multivalued velocity) arises in the classical formulation 
of the problem. To consider the asymptotic solution as the contact line is approached, we 
expand {p, u,v} = ^^${pi(9),Ui(9),Vi(9)} f l and find at leading order 

po{u + v ) = -p v , p =0, p = 0, 

subject to p' = on 9 = {0, ir}. The density is solved as po = 1/2, having imposed its 
expected value at the interface. To find the leading-order velocities, we continue to first order 
in the governing equations, where 

u = -v' , Ca k (vo + v' ) + pi + pi = 0, 
Ca fe (^ + t5 )-pi" -p[ = 0. 

with the wetting and no-slip conditions. We also assert that the profile must be flat at these 
very small distances to the contact line for a well-defined microscopic contact angle in the 
Young equation — requiring p(f,ir — 9s) = 1/2, at least up to this first-order correction, and 
leading to 

sin #s cos 9 + cos 9s sin 9 f u \ / — cos t 



Pi , . 

Vq I \ sm( 

We now consider the stresses and pressure, which in inner variables are scaled with Cn _1 
(readily seen from eqs. ©-([2])), as their singular behaviour in the classical model of |16] 
is the hallmark of the moving contact line problem. The total stress components in polar 
coordinates and in inner variables are 



Mff = Cafc 



'p\l-~pf + {d e p) 2 /f 2 -{d f p) 2 }/2 

+p{d 2 f p + d f ~p/f + d 2 e p/f 2 - p(l - p)(l - 2/5)}] 
+ p [4d f u/3 - 2 (5 + d e v) /(3f)] , 
M? e = -C&k^drp dgp/r + p [(d e u - v) jr + d f v] , 



M ee = Ca fc 



[p\l-~pf-{d e p) 2 /f 2 + (dr~p) 2 }/2 

+p{d?p + d f p/r + d 2 p/f 2 - p(l -p)(l~ 2/5)}] 
+ p [-2d f u/3 + 4 (u + d e v) /(3f)] , 

so by substituting in our results, we find at leading order M = O(l), as all 0(1/ f) terms 
cancel. To obtain the precise form of the stresses, the second-order terms in the governing 
equations are needed. The pressure in this inner region is given by p = p 2 (l — 3/5) (1 — 
p)/ (2 Cafe), so that as r —> it satisfies p = —1/(32 Cafe) + O(f), being finite at the contact 
line. 

Continuing with these second-order terms, we find the density and velocity corrections 
h = C h + Cp 2 cos(2<9), 



ui \ _ 

Vl 



sin0 s (cos(20) - l)/4 - C 5 sin(20) 
- sin 9 S sin(20) /4 + C s (l - cos (20)) 
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where the arbitrary constants Cp x , Cp 2 and C# would be set by the full solution of the inner 
problem. A possible flow scenario where Cp 1 = —0.1, Cp 2 = 0.3, C„ = — 1, for 9 s = vr/4 
is shown in Fig. [3ja). All of these results allow us to determine the leading-order stress 
components as 

(Mi cos(2#) - M 2 sin(20) + M 3 )/32 + 0(f), 
-(M 2 cos(20) + Mi sin(2#))/32 + 0(f), 
(M 2 sin(2#) - Mi cos(20) + M 3 )/32 + 0(f), 

where Mi = cos(20 s )/Ca fc +8sin0 s , M 2 = 32C e +sin(26> s )/Ca fe , and M 3 = (1 + 64(7^ ) /Ca fe - 
8sin#g/3, showing that the stresses are nonsingular as f — > 0. 

3 Extensions 

Having demonstrated the ability of our diffuse-interface model to alleviate the moving contact 
line problem with no-slip applied, we now consider a number of other features of contact line 
flow. A recent paper, [35], critically examined the interface formation model of Shikhmurzaev 
|34j . one of the more complex continuum theories proposed to deal with contact line flows. 
There, it was shown that the model degenerates to the same macroscopic flow as slip models 
but it was also seen to have features that most of the simpler models such as Navier-slip 
do not. In particular, the interface formation model captures finite pressure behaviour, it 
relaxes both the fluid-fluid and the fluid-solid interfaces through the modelling of surface 
layers, and the contact angle is able to vary dynamically from its static value (with its value 
determined rather than prescribed empirically). Finally, the fluid is able to 'roll', as in a 
moving frame of reference there is no stagnation point at the contact line, allowing particles 
to reach and transfer through the contact line in finite time. Whilst these features occur at 
lengthscales too small to probe with current experimental ability (see discussions in |35j), it 
is of interest that the diffuse-interface model is capable of similar predictions, with various 
features added. The model studied thus far already predicts finite pressure at the contact 
line, and alleviates the stagnation point predicted by slip models through mass transfer. It 
relaxes the fluid- fluid interface, but has classical no-slip at the wall in contrast to the effective 
slip of the interface formation model. Although not necessary, this may also be relaxed 
through carefully prescribing a generalised Navier boundary condition (GNBC), suggesting 
that the slip velocity is proportional to the total tangential stress (the sum of the viscous 
and uncompensated Young stress — arising from the deviation of the fluid-fluid interface from 
its static shape), and derivable using variational arguments from the principle of minimum 
energy dissipation [2U [25]. Our diffuse-interface model also prescribes the actual dynamic 
contact angle 9^ to be equal to the static value 9s through the wetting boundary condition. 
An alternative is for this condition to hold at equilibrium, with the density relaxing to it in 
finite time when out of equilibrium, as initially discussed (but not implemented) for binary 
fluids [17] , and more recently used in numerical simulations [24^ |2"5"[ l4"5[ l!"6] . 

For our liquid-gas configuration, the GNBC and generalised wetting boundary condition 
may be considered in an analogous fashion. In dimensional form the wetting boundary con- 
dition is generalised to 

a(dtp + u-Vp) = -L(p), (8) 



Mf f 
Mss 
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Figure 3: The asymptotic behaviour of density (contour plots) and velocity (arrows) as the 
contact line is approached. 

where L(p) = Kn w • Vp + is the wall chemical potential, and with a = representing 

instantaneous relaxation to equilibrium. The GNBC for this application with inverse slip 
length (3 is then 

L(p)(t w • Vp) - t nt = /3(u - u w ) • t w , 

with t w the tangent to the wall, and f n t the viscous shear stress. Note a = reduces to the 
popular Navier-slip condition. In nondimensional form in polar coordinates for steady flow 
in the inner region, eq. ([8]) reduces to 

L{p) = Ca fc II (u dfp + v i%/5/ f) , 

where L{p) = ±dep/r + p(l — p) cos 6s on 6 = {0, n}, and with a nondimensional parameter 
arising, II = ap^/^pie), describing the extent of the wall relaxation (II = being instanta- 
neous). Similarly the GNBC reduces to 

^L{p)dfp — Cafc p [dgu/f — v/r + dfv] = j3 Cafc(l ± u), 

where f3 = fie/pi is the nondimensional slip parameter, and along with II are both chosen to 
be formed with the interface thickness e such that they are considered as O(l) in the limit 
Cn^ 0. 

To consider how eq. © allows for microscopic contact angle variation dependent on flow 
conditions, we note that at equilibrium (denoted with subscript e) that eq. ([8]) gives L(p) = 
Kn w • Vp[ e + /^(pl/2) = 0. Now in a steady, dynamic situation (denoted with subscript 
d), eq. ([5]) implies a u • V/j| d = — Kn w • Vp| d — /4(pl/2), thus using the equilibrium result 
and considering this at the contact line with wall velocity u„ = — Vt w , we find VasinOa = 
—K(cos6d — cos 9s), or in nondimensional form 

Ca fc n = (cos 6*5 - cos^/sin^ « 9 d - 6 S , 



S 



where the final approximation holds for 9^ — 6$ <C 1, in agreement with the binary fluid case 

Another consideration is the behaviour near the contact line if density gradients near the 
wall are permitted far from the contact line. Our wall free energy was specifically chosen 
to prevent this, but we may also consider the two other situations used previously in the 
literature, namely (i) specifying a density at the wall p = p a on dQ, as in [23], and (ii) choosing 
a linear form in the density for the wall free energy f w (p) = ap, as in [311 HI 03] . To consider 
the contact line behaviour when (i) replaces the wetting boundary condition is straightforward, 
but for (ii), we must understand how to impose the microscopic contact angle to compare to 
our previous condition. Following [4], we use Young's law cos^s = (ay — aj^ja and compute 
ay and o~i, by integrating the free energy per unit area along the corresponding interface. 
This gives cos#s = [(1 — A) 3 / 2 — (1 + vl) 3//2 ]/2, where A = Aae/(Kpi) is nondimensional. 
This may then be inverted to give the appropriate value of A for a given contact angle 9s, 
and corresponds to the nondimensional boundary condition Cnn„, • Vp = —A/4. 

Adding these features into the diffuse-interface model do not dramatically alter the contact 
line behaviour, but subtle differences in the asymptotic results are demonstrated in Fig. [3] 
for selected arbitrary constants, where Ca^ = 0.1 and 9s = tt/4, and may be compared to 
the equilibrium situation in Fig. [2 The cases considered are (a) the original model without 
slip or wall relaxation, (b) using the linear form for the wall energy density in the wetting 
boundary condition, (c) adding a precursor film at the wall (where the wetting boundary 
condition is replaced by p = p a /pL = 0.53), (d) allowing finite wall relaxation, II = 5, (e) 
including the GNBC (J3 = 2) but with II = 0, and (f) the GNBC with = 2, U = 5. All 
models behave as expected, resolving the stress and pressure singularities, and including the 
effects they intend, e.g. capturing the film in (c), increased microscopic contact angle in (d) 
and (f), and reduced wall velocity in (e) and (f). There are only small differences between 
cubic and linear wall energy forms near the contact line, mainly that the linear form shows a 
broader band of density variation. This hints at the important difference that will occur near 
the wall but far away from the interface, where density gradients will remain present for this 
linear form. 

4 Conclusions 

We have shown analytically that a diffuse-interface model is able to resolve the moving contact 
line problem through relaxing the interface from being sharp to thin, without need to prescribe 
any additional physics at the contact line. Whilst slip, precursor films and finite-time wall 
relaxation have been considered, they are not necessary to resolve the moving contact line 
problem. We believe that the present study will motivate further analytical and numerical 
work with diffuse-interface models, such as to consider heterogeneous walls [281 EH E3 H2] . Of 
particular interest would also be the inclusion of non-local terms into the governing equations; 
this was considered in [21] for equilibrium wetting using a density-functional theory. The 
sharp-interface limit of binary fluids is also an interesting question. Preliminary investigation 
suggests that the models numerically considered in [251 [251 HZ] will exhibit the asymptotic 
problems encountered in [13] for general phase-field models without boundaries (as eluded to 
earlier), suggesting that the brief work of [52] should be considered only as a first attempt at 
this challenging problem. 
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